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ABSTRACT 

Outflows promote the escape of Lyman-a (Lya) photons from dusty interstellar media. The pro¬ 
cess of radiative transfer through interstellar outflows is often modelled by a spherically symmetric, 
geometrically thin shell of neutral gas that scatters photons emitted by a central Lya source. Despite 
its simplified geometry, this ‘shell model’ has been surprisingly successful at reproducing observed 
Lya line shapes. In this paper we perform automated line fitting on a set of noisy simulated shell 
model spectra, in order to determine whether degeneracies exist between the different shell model pa¬ 
rameters. While there are some significant degeneracies, we find that most parameters are accurately 
recovered, especially the HI column density (Am) and outflow velocity (u exp ). This work represents 
an important first step in determining how the shell model parameters relate to the actual physical 
properties of Lya sources. To aid further exploration of the parameter space, we have made our 
simulated model spectra available through an interactive online tool. 

Subject headings: radiative transfer - ISM: clouds - galaxies: ISM line: formation - scattering - 
galaxies: high-redshift 


1. INTRODUCTION 

Lyman-a (Lya) emission enables us to both find and 
identify galaxies out to the highest redshifts. Whether 
or not a distant galaxy has a Lya emission line appears 
to be clo sely connected with the presence or absence of 
outflows ( Kunth et al.|1998 Atek et aL|2008 ). Partially 
coherent scattering ot Lya photons oft the outflow pro¬ 
vides an efficient way of shifting them into the wing of 
the line profile, where the galaxy is optically thin. Re¬ 
cent near-IR spectroscopic measurements have confirmed 
that the Lya spectral line is systematically redshifted 
with respect to other non-resonant nebular lines 
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the Lya spectra emerging from galaxies. The simulations 
can broadly be split into two categories: 
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Smith et al. 2015). It is worth stressing that this 
is a challenging task as simulating the distribution 
and kinematics of neutral gas in feedback driven 
outflows - which play key roles in the Lya transfer 
process - require a sub -pc spatial resoluti on (com¬ 
pare e.g. figs. 10, 11 in|Fuji taet al.| [2009). 


• The second approach is to simplify the complex 
topography of galaxies enormously, and use an ab¬ 
stract geometrical setup instead. This approach 
is computationally cheap, which allows model Lya 
spectra to be generated for a large set of parame¬ 
ters. 


A particular example from the second group, which has 
been very successful i n reproduc i ng observed Lya spec¬ 
tra, is the shell mode l ( Ahn|2004 Verhamme et al.||2006 
Schaerer et al. 20111. This model consists of a centra 


Lya source surrounded by an outflowing thin, spheri¬ 
cal shell of hydrogen and dust. The shell model has six 
model parameters: the equivalent width of the emitting 
source, EWf, the intrinsic width of the Lya spectrum, 
cr,; the hydrogen column density, JVhi; the dust optical 
depth, T(i ; the temperature, T; and the outflow velocity, 

^exp* 

Due to its simplicity, and the ability to reproduce a 
number of observed spectra, the parameters of the shell 
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use a discrete set of models from which they choose the 
best fitting spectrum based on y 2 values. A systematic 
study of shell models and their physical relevance has yet 
to be performed, however. A number of open questions 
include the following. 
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1. Is there a unique mapping between shell model pa¬ 
rameters and Lya spectra, or can several different 
parameter combinations - with different physical 
implications - produce indistinguishable spectra? 
If the latter is true, two equally good fits to an 
observed spectrum could be found that have very 
different physical meanings, and so caution must 
be exercised when drawing conclusions from shell 
model parameter fits. 

2. How uncertain are the parameters obtained from 
the model fitting procedure? The best-fitting shell 
model parameters to observed spectra are often 
quoted without estimates of uncertainties. 

The two questions are clearly related to one another: 
question (1) focuses on the intrinsic degeneracies between 
parameters, while question (2) concerns specific example 
spectra with some measurement uncertainty, and the ef¬ 
fect this has on the accuracy with which the model pa¬ 
rameters can be recovered. Both must be understood if 
we are to use the observed spectral line shapes of Lya 
to infer the physical properties of outflowing gas. This 
is especially important if one considers that, at high-z, 
Lya is likely the only line that will allow us to put any 
constraints on outflows of atomic hydrogen gas. 

The goal of this paper is to address these questions, and 
to develop a procedure to automatically fit shell model 
parameters to ‘realistic’ (i.e. noisy) Lya spectra. The 
paper is structured as follows. In Sec. [2] we explain our 
method for constructing simulated Lya spectra, and de¬ 
scribe the fitting procedure we use. We present the re¬ 
sults of applying this procedure to a range of simulated 
spectra in Sec. [3j and discuss the implications of our re¬ 
sults in Sec. [4j We then conclude in Sec. [5] 

2. METHOD 

We brief ly d escribe the Lya radiative transfer calcu¬ 
lations in §2.1[ and present our application to the shell 
mod el in f j2.2| Our fitting procedure is then described in 

In the following, log denotes the logarithm to base 10, 
In is used for the natural logarithm. The wavelength 
is given in units of velocity offset from line center, v = 
c(A/Ao — 1), where Ao ~ 1215.67A corresponds to the 
wavelength at line centre. 


2.1. Lya radiative transfer 

We use the Lya Monte Carlo (MC) ra diative transfer 
code tlac. w hich was previously used in |Gronke fe Di-| 
jkstra (2014). We briefly summarize the main points of 
the algorithm here. For a m ore detailed review we refer 
the reader to Dijkstra (2014). 

In a radiative transfer MC code, the photons are repre¬ 
sented by individual particles (or rather, photon packets) 
whose paths of propagation in real and frequency space 
are tracked throughout the simulation domain. We re¬ 
peat the following steps for each photon in the Monte- 
Carlo simulation: 


1. A photon is emitted in direction fc;, with intrinsic 
frequency Vi drawn from some Lya source emission 
pdf, f(v). 


TABLE 1 

Simulation grid for shell model parameters 


Symbol 

Values 


Units 

Vexp (a) 

(0, 2, 5, 8,10,15, 20, 30,40,..., 490) 

kms" 1 

log Aim (b) 

(17.0,17.2,... 

. ,21.8) 

log(cm —2 ) 

logT W 

(3.0, 3.4, 3.8.. 

..,5.6) 

l°g(K) 


[1, 800] 

(continuous) 

kms -1 

Td (6) 

[0, 5] 

(continuous) 

- 

EWi (f) 

[1, 150] 

(continuous) 

A 


Cl expansion velocity, Cl hydrogen column density, C) (effective) 
temperature, widtlri of intrinsic spectrum, Cl dust optical depth, 
Cl intrinsic equivalent width. 


2. A number r is drawn from an exponential distri¬ 
bution, quantifying the distance d the photon will 

travel. Here, d is given by r = /(^(ohiUhKs) + 
Td«d(s))ds, where nx and ax represent the num¬ 
ber density and the cross section of dust ( d ) or hy¬ 
drogen (hi), respective!^ 


3. The position of the photon is updated to p —> 
p + kd. At the new position, the photon 
is either scattered by hydrogen (probability of 
nmo-Hi/(o-Hi«Hi(s) + cr d nd)), or absorbed by dus!0 
In the former case, a new direction is drawn from 
the proper phase function (i.e. the angular redis¬ 
tribution function) and a new frequency is drawn 
from the frequency redistributi on function (which 
depends on direction, see Dijkstra & Kramer|2012). 


4. Steps (2) and (3) are repeated until the photon 
escapes the simulation domain or is absorbed. 


If the photon escapes the simulation domain, its fre¬ 
quency and other properties are recorded. This proce¬ 
dure is repeated until the desired total number of pho¬ 
tons is reached. Their combined frequencies then yield 
the simulated Lya spectrum. 


2.2. Implementation of the shell model 

As discussed in Sec.[T] the shell model parameters con¬ 
sist of: (a) the source properties - namely, the equivalent 
width of the emitting source EWi, and the intrinsic width 
of the Lya spectrum erg and (b) the shell properties, i.e. 
its hydrogen column density Ahi, its dust optical depth 
Td along a path that passes radially through the shell (the 
dust is confined to the neutral shell), temperature T, and 
the outflow velocity u eX p- In order to cover a particular 
domain of the parameter space with as few simulation 
runs as possible, we split these six parameters into two 
sets: the simulation parameters (u exp , Am, T), and the 
post-processing parameters ( EWi , cq, Td)- The simula¬ 
tion parameters change the spectrum in a complex way. 
We therefore discretize these parameters on a grid and 

1 Note that am depends strongly on frequency, which generally 
translates to a dependence of am on s. In additio n, a hi is a 
function of the temperature T (see, e.g., Dijkstra|2014| for details). 

2 For simplicity, we assume all-absorbmg (i.e. non-scattering) 
dust throughout this work (that is, we set the dust albedo to 
A = 0). We therefore have r a = (1 — A) r,/ = 77 /, where r a 
[r,j] denotes the optical depth to dust absorption [absorption & 
scattering]. Thi s as sumption does not affect our predicted spec¬ 
tra (Laursen et alT)|2009|, provided that results are ‘rescaled’ as 

—> t ( ] / ( I — A) when comparing to models with non-zero A. 
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run our radiative transfer simulation at each combina¬ 
tion of parameter values. This results in agrid of Lya 
spectra with parameters presented in Table [T| The exact 
discretization was chosen so that the grid is fine enough 
to allow for an accurate likelihood analysis, i.e. so that 
the likelihood functions of our simulated spectra are suf¬ 
ficiently well sampled. 

The post-processing parameters do not require explicit 
simulations, but can be modified a posteriori. This con¬ 
siderably reduces computational cost, and allows the pa¬ 
rameters to be varied continuously (i.e. they are not 
restricted to discrete values on a grid). Each photon has 
its own unique weight in our simulations which is given 

by 


_ fp( y i) - T d N H 1 /N HI 
fr(Vi) 


( 1 ) 


where Vi is the velocity at which we inject the photon 
and .ZVhi is the hydrogen column density actually encoun¬ 
tered by the photon in the simulation”] Furthermore, f r 
denotes the emission pdf that we used in our Monte- 
Carlo simulation (characterized by a Gaussian emission 
line with a = 800 kms -1 plus continuum, see below), 
while f p denotes the emission pdf that we wish to sim¬ 
ulate in post-processing (which is also a Gaussian with 
standard deviation er,; plus continuum). This means f r 
and f p are of the form 

/(».) = {fd£ + s ( >'<M f«M<¥ (2) 

I 0 otherwise 


where G(v\a) denotes a normalized Gaussian with stan¬ 
dard deviation a centered at u = 0, An is the bandwidth 
of the spectrum, and = AuAl ya /(cEWi), where AL ya 
is the emission wavelength of Lya. We stress that each 
photon is assigned an individual weight. This means, the 
post-processing parameters also affect the shape (and not 
only the normalization) of the resulting spectrum. We 
tested our method by comparing several post-processed 
spectra with ‘real’ spectra where Td, a l and EWi were 
given as fixed input parameters in our Monte Carlo ra¬ 
diative transfer code. 

Note that in this framework the escape fraction of Lya 
photons /e S c can be calculated using 


1 

N 


E 

photons 


G{ViW,) -tN hi /N hi 
fr{Vi) 


( 3 ) 


where Af = YjG{ v iWi)/f r (vi). 


and observational studies (e.g. Blanc et al. 

2011 

Hayes 

et al. 

2011 

Dijkstra & Jeeson-Damel 

2013 

and can be 


measurements are available, see (4.1). Alternatively, the 
obtained f esc (and its uncertainties) can be seen as an 
additional result of the shell-model fitting. 

In order to minimize the uncertainty on the post- 
processed spectrum, the ratio f p {vi)/f r (vi) should be as 
close to unity as possible across all frequency bins. When 
performing the grid simulations, we therefore chose a 


3 This column density JVjji > IVhi because scattering increases 
the total path a Lya photons traverses before it escapes from the 
shell. 


‘typical’ fiducial emission pdf, f r (vi ), with fixed values 
of (uj, EWi) = (800 kms -1 , 2.92A) for all simulations. 
The large cq and small EWi ensure that a wide range of 
v bins is initially populated with photons in the simula¬ 
tions, so f r > 0 in any birr] The simulated value of dust 
optical depth is Td = 0 lor minimizing computational 
time. 

In our analysis, we allow the post-processing param¬ 
eters to vary continuously between (cq, Td, EWi ) = 
([1, 800] kms -1 , [0,5], [1, 150]A). These ranges were 
chosen to try and capture all physically-interesting sce¬ 
narios, but are still somewhat arbitrary. They can be 
extended if necessary without requiring additional com¬ 
putations. 

The resulting dataset consists of 10,800 spectra over a 
grid of the 3 simulated parameters, with 170,000 photons 
per spectrum. The number of photons was chosen so 
that the approximate Poisson uncertainty (for a uniform 
spectrum) is < 5% on each bin; \/iVbins/170, 000 ~ 2.4% 
(-Wins = 100). The spectra can be accessed through an 
interactive online too@ 

2.3. Simulated spectra and fitting procedure 

To test the effectiveness of spectrum fitting in re¬ 
covering shell model parameters, we produced a fully- 
simulated dataset consisting of 50 spectra for a range 
of randomly-chosen input parameters]^] The simulated 
spectra w ere f ound directly by using the method de¬ 
scribed in §2.1[ i.e., no post-processing was used, yielding 
the ‘true’ (input) spectrum, F(vi;0 0 ), for each set of in¬ 
put parameters 9q, and then adding noise, rq, to each 
spectral bin i, resulting in an ‘observed’ spectrum 

di = F{vi, e 0 ) + m. ( 4 ) 

While in reality the noise level in each bin will be a com¬ 
plicated function of instrumental characteristics, observ¬ 
ing conditions, and the properties of the data reduction 
pipeline, for simplicity (and to avoid tying our results 
to a particular observational configuration) we chose n* 
to be Gaussian, with zero mean, constant variance, and 
no correlations between frequency bins. The frequency 
binning was chosen as constant intervals of width A= 
50 kms -1 over the rang^Ju = [—2500, 2500] kms -1 , with 
a simple tophat bandpass and no smoothing. This cor¬ 
responds to a relatively high (but plausible) resolving 
power of R = c/Avi ~ 6000, which lies within the ca¬ 
pabilities of, e.g., the IRCS instrument on the Subaru 
telescop ep| We discuss the impact of changing the bin¬ 
ning in §T2| 

4 A uniform emission pdf would achieve the same effect, but this 
would be inefficient, overpopulating the wings of the spectrum. 

5 The online tool can be accessed under http://maxgronke. 
wordpress.com/tools/tlac_web/ or 

http://bit.ly/man-alpha On the website, it is possible to (i) 
plot up to tour shell-model spectra, (ii) upload and display own 
spectra, and, (Hi) download the plotted data. 

6 These were drawn from a uniform distribution with the bounds 
given in Table [T) except for r^, which was drawn from [ 0 , 2 ] to save 
on computational resources. 

7 While we expect this range to be sufficient for most observed 
spectra (i.e. wide enough to sample the L ya line and the adjacent 
continuum), it is possible to alter this range if needed since our 
method is not dependent on this choice. 

8 http://subarutelescope.org/Observing/Instruments 
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Fig. 1.— Examples of similar spectra produced by significantly different sets of shell model parameters. The ‘observed’ spectrum is shown 
in gray, the spectrum obtained from the true shell-model parameters in black, and one of the best-fit models produced from 640 trials of 
the fitting procedure is shown in blue. The spectra were chosen to highlight extreme examples, one for e ach o f the six model parameters 
(highlighted in bold in each panel). In each case, the fitted model is a better fit than the true model; see j |3.1| for details. The parameters 
have the following units: (^exp, -/Vhi, T, a*, EWi) = (kms -1 , cm -2 , K, kms -1 , A). 


Following our assumption of uncorrelated Gaussian 
noise, the appropriate goodness-of-fit statistic for the 
spectra is 



F{vj\9) 

Vi 


2 


(5) 


where the sum is taken over the frequency bins. This de¬ 
fines the likelihood C(0\d) oc e -Ax / 2 . We systematically 
explore the likelihood surfaces of our simulated data us¬ 
ing Markov Chain Monte Carlo (MCMC) and nonlinear 
optimization (fitting) methods in the next section. Since 
we have assumed that the noise rms (i.e., the statistical 
weight) <7 i = &(vi) = const., and is independent of model 
parameters 6, the statistic 


A 2 = Y,(di-F(vi-,e)) 2 ( 6 ) 


can be used interchangeably with Ay 2 by optimization 
algorithms, a fact that we will also use in Sec. [3] 

Note that there is also a ‘theoretical uncertainty’ asso¬ 
ciated with the model spectrum F(v;0). Because each 
spectrum is calculated through a Monte Carlo radiative 
transfer procedure with a finite number of photons, there 
is a Poisson error on the value of F in each bin, deter¬ 


mined by the number of photons that ended up in that 
bin (and their respective weights). Since the number of 
photons per bin also varies as a function of the input pa¬ 
rameters, the theoretical uncertainty will vary across the 
parameter space in a potentially complicated way. This 
uncertainty has been mitigated as much as possible by 
our use of a large number of photons to calculate each 
model on the grid of simulations, as described in §2.2| 
This helps to ensure that the assumed observational un¬ 
certainty, (jj, dominates the Poisson uncertainty in each 
bin. Under these conditions, the theoretical uncertainty 
can safely be ignored, and Eq. © remains valid. 


3. RESULTS 

3.1. Intrinsic degeneracies 

To investigate whether the best-fit parameters recov¬ 
ered by the standard spectral fitting procedure could be 
misleading, we performed 640 local A 2 -minimizations on 
each of the 50 spec tra in our fully-simulated observa¬ 
tional dataset (see § 2.2) without adding any additional 
noise (n* = 0). We chose the number of local minimiza¬ 
tions rather arbitrarily, but sufficed to show a few spe¬ 
cific cases that highlight the dangers of blindly following 
this procedure. We started each minimization from a 
randomly chosen initial guess. We then discarded the 
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Fig. 2.— Example analysis of an asymmetric double-peaked profile. Upper right : Simulated data (black solid line) with its associated 
68% CL observational uncertainty (gray shaded region), and the spectra for the true (input) and median MCMC-estimated parameters 
(light red and light blue, respectively). Main plot : The ID histograms and 2D contour plots show the one- and two-dimensional marginal 
posterior distributions from the MCMC chains respectively. The red solid line and red marker show the true (input) model parameters. 
In the histograms, the dashed lines mark the 16th, 50th and 84th percentiles; these numbers are also stated on top of each column. In 
the contour plots, the blue lines mark the (2, 1.5, 1, 0.5) cr contours, and the gray shading gives the posterior density. Recall that the 
parameters log A/hi, logT and u eX p are discrete (visible through the blocky structure of the contours). 


fits where A 2 was less than the A 2 for the true model 
parameters. A more systematic study t hat i ncludes ob¬ 
servational uncertainties is presented in §3.2| 

Fig.Hshows examples of degenerate model fits for each 
of the shell model parameters. In each panel we display 
the ‘observed’ spectrum (in gray), the spectrum obtained 
for the true input model parameters (black), and an al¬ 
ternative best-fit model (blue). It can be seen that no 
shell model parameter is protected from catastrophic re¬ 
covery errors, as the magnitude of the difference in each 
of the example cases allows significantly different physi¬ 
cal interpretations. 

As a concrete example, the central panel in the bot¬ 
tom row of Fig. [I] shows a spectrum that can be ex¬ 
plained with a very high dust content (t^ = 5), although 


it was actually produced by a model with very little dust 
{Td = 0.29). In fact, the value of = 5 here is the (ar¬ 
bitrarily cbosen) upper limit of our parameter space (see 
|2.2[ ), so an even higher value may be possible. In this 
model, the low value of cr,;, and the high outflow velocity 
combined with the low hydrogen column density, results 
in essentially no interaction between the photons and the 
hydrogen. The path length through the shell (and thus 
the chance of absorption) is therefore approximately the 
same for each photon, and so has little or no effect 
on the shape of the resulting spectrunfi] The same phe- 

9 In the notation of Eq. §• this means that in this particular 
case TVhi ~ Nhi for all photons leading to no spectral shape change 
due to a change of t^. 
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t) ex p (kms x ) = 9.65^7 


— “Data” 

- True 

Median 



* V ® ® ‘b > ^ ^ 

«exp (kms -1 ) log JVki/cm -2 logT/K a,-(kms -1 ) 


# <$> ^ 
EW, (A) 


Fig. 3.— Same as Fig. [2] but for a simulated spectrum with a Lya absorption feature. 


nomenon also explains why a ~ 1.5 order of magnitude 
decrease in HI column density barely affects the spec¬ 
trum in the central panel of the top row. 

Another case where a parameter hits its prior bound¬ 
ary is in the lower left panel of Fig. [l] where <r* reaches 
its maximum permitted value of 800kms -1 . The reason 
the intrinsic spectral width plays no role in the shape of 
this absorption feature is because the high hydrogen con¬ 
tent causes each Lya photon to be scattered many times. 
The relatively large value of r d leads to a high absorp¬ 
tion probability, yielding this broad absorption feature. 
While in this case <jj plays only a minor role, with dif¬ 
ferent values leading to qualitatively similar spectra, the 
same goodness of fit can only be reached by also changing 
the other parameters, in particular r d and EWi. This can 
be explained as follows: while the increase in cq means 
that there are initially more Lya photons in the wing, 
which have a higher chance of not being absorbed, this 
is compensated somewhat by the higher dust content. 


Since this increase also (and more strongly) affects pho¬ 
tons travelling longer distances (i.e. the ones closer to the 
core), a higher number of Lya photons is also needed to 
return a similarly good fit as the true model. 

The EWi — T~d degeneracy is also observable in the 
bottom-right panel, although the difference in EW; is 
not as large. This is because, for emission features, EW, 
seems to be reasonably well constrained by the contin¬ 
uum level. 

Finally, the expansion velocity is thought to be well 
constrained by the position of the spectral peaks and 
troughs. The situation shown in the top-left panel is 
therefore less straightforward the effects of changing of 
v exp can only be compensated by jointly changing EWi , 
tJi, and T d . 

3.2. Parameter uncertainty 

We now turn our attention to the uncertainty asso¬ 
ciated with the estimated parameters obtained through 
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Fig. 4.— Recovered (y axis) vs. input (x axis) shell model parameters. The error bars mark the 68% CL, between the 16th and 84th 
percentiles. The colors and symbols indicate different types of spectra: absorption features (red discs), P Cygni-like profiles (blue squares), 
single-peak (green triangles), and double-peaked (purple diamonds). 


shell model fitting. First, we consider two example cases, 
each represented by a spectrum for which the input shell 
model parameters were correctly recovered. Note that, 
our simulated dataset contains absorption, single- and 


spect ra 


a (e.g. Steidel et al. 

2010 Rivera-Thorsen et al. 

Yang et al. 2015). 

(Jut of these, the two cases 


absorption profile to show the extent of possible degen¬ 
eracies. Moreover, in some cases (especially at higher 
redshift) the blue peak may be further suppressed by 


the intergalactic medium (Dijkstra et al. 
al. 


et al. 112011 1, in which case our double peaked example is 


2007 [L 


jaursen 


actually more representative of what has been observed. 

We begin by constructing a si mul ated observation of 
each spectrum, as described in §2.3| The noise in each 
bin is chosen to have an rms of a~i = 0.5/, where / is 
the mean intensity of the spectrum. Note that while the 
noise properties have been chosen rather arbitrarily, the 
procedure does not depend on this choice, and also works 
with more realistic data errors. This choice of ch ensures 


that ai dyPoisson (see 12.31, and can be thought of as 
representing random errors due to instrumental effects 
and so on. 

This noise level corresponds to signal-to-noise ratios 
(SNRs) of ~ 15 — 50, which is c omparable with existing 
surveys (e.g. Adams et al. |2011[ ). We estimated the SNR 
by maximizing the quantity JT of, where the 

sums are taken over several adjacent bins. This corre¬ 
sponds to the standard procedure in observing pipelines 
(S. Wilkins, priv. comm.). For double-peaked profiles 


and absorption features, different measures are some¬ 
times used - e.g. taking the difference between the con¬ 
tinuum level and the absorption feature - that tend to 
result in higher SNR estimates. We will not consider 
such measures here. 

Next, we sample the likelihood of Eq. ([ 5 ]) using the 
affine invariant Markov chai n M o nte Carlo (MCMC) 


ensemble sampler emcee ( 

Goodman & Weare 

2010 

Foreman-Mackey et al. 

2012 

). We use 900 walkers with 


starting p osit ions of the walkers, we used the A-minima 
found in §3.1[ weighted by the value of the likelihood at 
that point, plus a small random perturbation to avoid 
producing initial paths that are too similar. 

Fig. H shows the results of the MCMC parameter 


estimation for the first example casq__| an asymmet¬ 
ric, double-peaked profile with an estimated SNR of 
~ 32. The true (input) shell model parameters are well- 
recovered, falling within the 68% credible interval for all 
but Td, which is estimated to be slightly higher than its 
actual value. One can also see that the expansion ve¬ 
locity and column density are very well constrained (in 
fact, reaching our grid resolution limit), and that the lcr 
uncertainty on the temperature is almost half an order 
of magnitude. 

The second example, a spectrum with a broad ab¬ 
sorption feature and SNR ~ 20, is shown in Fig. [3] The 
uncertainties on the estimated shell model parameters 

10 This and the other triangle plots were produced using a mod- 


“ this and the other trian gle pl ots were produced us ing ; 
ified version of triangle.py (Foreman-Mackey et al. 2014). 


















































































































are much larger than in the previous example, with all 
but A'hi being relatively poorly constrained. Degenera¬ 
cies between many of the parameters can be seen clearly; 
the values of a,;, EWi, A'hi and are all correl ated 
with one another, for the reasons discussed in §3.1[ 
Note that the distinct directions visible in the — EWi 
plane, and other ‘striping’ features in the plots, are a 
result of the discrete parametrization of log A'hi , which 
is mostly localized to the three preferred grid values of 
[20.4, 20.6, 20.8]. This was confirmed to be the case by 
plotting the MCMC steps falling in these three column 
density bins individually - it is not an artefact of the 
sampling procedure. A higher-resolution gridding of the 
log A'hi parameter would reduce this effect. 


Fig. a shows an overview of the estimated parameters 
for all 50 spectra in our simulated dataset, after run¬ 
ning exactly the same MCMC procedure as for the two 
example spectra. The data points are colored accord¬ 
ing to spectral type, which was assigned to each spec¬ 
trum following visual inspection. The parameter esti¬ 
mation procedure reliably recovers at least some of the 
shell model parameters for all spectral types, particu¬ 
larly the expansion velocity u exp , and column density 
A'hi > 10 19 cm -2 (although the uncertainty is larger for 
smaller column densities). The intrinsic source parame¬ 
ters cTj and EWi are also mostly well-recovered, but with 
a bigger uncertainty. An exception is for spectra with 
absorption features, where ct,; is typically overestimated 
by ~ 200 kins -1 . The worst constraints are obtained for 
T and r^, for which the scatter and uncertainty often 
reach the prior boundaries for those parameters. 

4. DISCUSSION 


Given the large parameter space spanned by the shell 
model parameters, and the consequent variety of Lya 


spect ra (as is nicely illustrated in fig. 5 of Schaerer et al. 


2011), it is difficult to make foolproof statements about 


parameter degeneracies and uncertainties. We generally 
recommend comparing individual spectra to the theoret¬ 
ical predictions on a case-by-case basis instead. In this 
Section we discuss how we can nevertheless reach more 
generic conclusions on parameter degeneracies and un¬ 
certainties by analysing a sample of 50 mock spectra. 


4.1. Parameter uncertainties and degeneracies 

From the results of our likelihood analysis, we find 
that: 


• Lya spectra are very sensitive to the expansion ve¬ 
locity, v exp . In particular, double-peaked profiles 
allow the underlying velocity field to be recovered 
with great accuracy (limited by the simulation grid 
spacing to ~10kms -1 in this study). Absorption 
features and single peaks lead to greater uncertain¬ 
ties of up to ~ 50kms -1 . This is con sistent with 
the findings of Verhamme et al. (2008). 


• The intervening column density, A'hi, can also be 
recovered well - at least for A'hi > 10 19 cm -2 . In 
this high column density regime, the (68% CL) 
uncertainty was limited by our grid spacing to 
A log A'hi /cm -2 = 0.2, while for lower column den¬ 
sities the uncertainty can be up to one order of 
magnitude. 


• The effective temperature (including turbulent mo¬ 
tion as well as the true temperature) has a complex 
effect on the Lya spectrum, and in most cases this 
parameter cannot be usefully constrained (also see 
e.g. Verhamme et al. 2008). 

• The width of the intrinsic spectrum, cr,, can be re¬ 
covered for Lya emission features with very high 
accuracy (Act,; < 50kms -1 ). For absorption fea¬ 
tures, on the other hand, this parameter is not con¬ 
strained so well (A Vi ~ 200kms -1 ). 

• We found that the dust content, t>i, was essentially 
unconstrained by the shape of the Lya spectrum. 
Exceptions to this are for some double-peaked fea¬ 
tures with low dust content, where the dust content 
defines the ratio of the heights of the peaks. If, 
however, measurements of the Lya escape fraction 
/e S c are available, it is possible to use this extra 
information in our analysis. When then expect 

to be better constrained in some cases. 


• The intrinsic equivalent width, EW, , can mostly be 
recovered with an uncertainty of ~ 20A regardless 
of spectral shape and the value of EW t . 

Note that these numbers depend on the spectral resolu¬ 
tion and quality (SNR). We present the result of changing 
these quantities in the next section. 

We also want to highlight that additional parameters 
can easily be included into our analysis which might lead 
to further degeneracies. For example, the line-center of 
the intrinsic spectrum can be unknown or afflicted with 
uncertainties when dealing with real data. We found that 
in some cases this ‘spectral shift’ is degenerate with the 
outflow velocity n eX p and/or the column density A'hi- 


4.2. Impact of SNR and frequency binning 

To test the sensitivity of our results to the assumed 
measurement uncertainties, we randomly selected 10 
out of th e 50 simulated spectra and repeated the anal¬ 
ysis of [3.2 for five different noise levels, Si/I = 
(0.25, 0.5, 1, 2 4) We also include the two spectra 
analyzed in §3.2[ for which the resulting SNRs are 
~ (64, 32, 17, 9, 4) and (39, 20, 10, 6, 3) respectively. 
Note that different noise realizations were drawn for each 
of the noise levels, for each spectrum. 

Fig. m shows the resulting parameter uncertainties as 
a function of SNR for each of the 10 + 2 spectra. For 
Vexp , A'hi, and cr,;, there is a flattening in the uncertainty 
beyond SNR ~ 20 for all types of spectrum. A similar 
behaviour is found for EWi but at a slightly higher SNR 
of ~ 30, although in this case the curves do not become 
completely flat, even at the lowest noise levels we con¬ 
sidered. For the remaining two parameters (r^, T), the 
evolution of the uncertainty with increasing SNR is not so 
systematic, with different behaviors for different spectral 
types. For T, some of the spectra show the same flatten¬ 
ing at high SNR as the other parameters, for example, 
while others are considerably more scattered. Note that, 
because the noise realizations change between the differ¬ 
ent choices of noise level, it is not surprising that there 
should be some random variation in the goodness of fit 
and parameter uncertainties. For there is a split be¬ 
tween spectra for which the uncertainty rapidly decreases 
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Fig. 5.— Parameter uncertainties vs. signal-to-noise ratio (SNR) for 12 different mock spectra, with five choices of d, for each. The y 
axis shows the 68% CL error estimate for each parameter, calc ulate d as the difference between the 84th and 16th percentiles of the sampled 
posterior distribution. The SNR is calculated as described in §3.2| 


with SNR, and those for which it changes only slightly 
(presumably in cases where the spectra are relatively in¬ 
dependent of Td, so it remains poorly constrained). 

In summary, for SNR ~ 30, four of the six shell 
model parameters can be obtained up to an un¬ 
certainty of (At) exp , A log -/Vni/cm -2 , Acq, A EWf) ~ 
(70kms _1 , 0.5, lOOkms -1 , 40A) for the majority of 
spectra, and these limits can be further decreased by 
increasing the SNR. 

Finally, we also varied the frequency binning (A m = 
(25, 50, 100, 200, 250, 300, 500) kins” 1 ) for several of 
the mock spectra, in order to test the influence of the 
spectral resolution on the parameter estimates. Fig- 
ure[6]in the Appendix shows that our constraints depend 
weakly on R provided that R >10 3 . For lower R the con¬ 
straints on several parameters (especially v exp ) increases 
rapidly for R < 10 3 , but only for some models. For ex¬ 
ample, for the case presented in Fig. [2] the characteristic 
double peak disappears for Avi = 250kms _1 (R ~ 10 3 ), 
which naturally leads to drastically increased uncertain¬ 
ties on many of the parameters. We therefore caution 
against too general conclusions on the precise impact of 
R , as this strongly depends on the spectrum that is an¬ 
alyzed. 


outflows, then the shell model provides an extremely 
quick and useful way to capture information about 
radiative transfer processes on scales that are difficult 
to model from first principles. 


As a first step towards addressing this question, we 
have constructed a large library of shell model spectra 
and which are available through an interactive online 


tool. A similar library was presented by Sc haerer et al. 
d201 1|) which was used in fitting observed Lya spectra 
(e.g., Hashimoto et al. 20151. Our work differs in that our 
grid of spectra is sufficiently closely-meshecpl to allow a 
systematic study of possible degeneracies between model 
parameters, via a fully automated likelihood analysis. In 
order to make this computationally tractable, we used 
a post-processing technique to model three of the six 
shell model parameters. We then applied an automated 
fitting routine to a small simulated dataset of ‘noisy’ 
shell model spectra, to see if a likelihood analysis can 
accurately recover the true (input) model parameters. 
This presents an important test of the reliability of 
shell model fitting procedures that are being used in 
the literature, and provides a simple estimate of the 
uncertainties on the recovered best-fit model parameters. 


5. CONCLUSION 

The ‘shell model’ for Lya transfer through galactic 
outflows has enabled a wide range of Lya spectra to be 
characterised in terms of six parameters. An important 
goal is to understand what physical information is actu¬ 
ally encoded in these parameters. This is an important 
question: if the parameters of the shell model contain 
a connection to the underlying physical properties of 


Our main findings include that from the spectral shape 
alone it was not possible to provide good constraints 
on the dust content, 7q, and effective temperature, T, 

11 Quantitatively, we have s ampled the funy,,. log Nvn, log T) 
space with 10800 models, while Schaerer et al. (2011) have 780 
models to cover this parameter space. In addition, we a llowed the 
dust opacity to vary continuously while |Schaerer et al. ([2011) ran 
8 discrete dust opacities. 
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in most models. In contrast, we were able to recover 
especially the outflow velocity (u eX p) and the hydro¬ 
gen column density (IVhi) well (especially when IVhi > 
10 19 cm -2 ). These two results reflect the fact that Lya 
radiative transfer observables are most sensitive to the 
hydrogen distribution and kinematics. Dust affects these 
observables (emerging spectrum, total flux), but how 
much depends on other parameters such as AVii , u exp 
and cq. This is consistent with observations which indi¬ 
cate that Lya escape is greatly affected by gas kinemat¬ 
ics, and which is likely responsible for the observed large 
scatter between f psr and Td (see, e.g ., Kornei et al.|201lf 
Blanc et al.||2011 Hayes et al.||2011 ). 


More quantitatively, lor spectra with signal-to- 
noise ratios SNR ~ 30 and spectral resolving 

power R = 6000, four of the six shell model 

parameters can be obtained up to an uncer¬ 
tainty of (Auexp, AlogiVni/cm -2 , Acq, AEWi) ~ 
(70kms _1 , 0.5, lOOkrns -1 , 40 A) for the majority of 
our spectra (see Fig [5] & § |4 1| ). These constraints 
can be improved by increasing the SNR (Fig [5]). We 
found a weak dependence of our results on R down 
to R ~ 1000 (below which our constraints deteriorate 
faster). However, statements about the impact of the 
resolving power R depends in detail on the features in 
the spectrum (e.g. the presence of double peaks clearly 
can depend on R). 

In general, precise quantitative statements on the 


magnitudes uncertainties and/or the extent of de¬ 
generacies depends on the location in the parameter 
space. (For example, it was not possible to recover eq 
specifically for spectra with an absorption feature.) We 
therefore recommend a case-by-case analysis of observed 
spectra using a f ull l ikelihood analysis (similar to the 
one presented in (3.2|. 


The likelihood analysis allows us to characterise spec¬ 
tra in a systematic, automated way, and produce robust 
estimates of parameters, their uncertainty, and the de¬ 
generacies with other spectra. This procedure will be 
useful in the near future, when a larger sample of Lya 
spectra with sufficient SNR and spectral resolution be¬ 
comes available. In a follow-up paper (Gronke et al. in 
prep) we will generate spectra for more realistic models, 
and fit shell models to these. This analysis will help us 
take the next & final step towards understanding what 
physical information is contained in the shell model pa¬ 
rameters. 


We thank the anonymous reviewer for useful com¬ 
ments. M.G. thanks Llufs Mas-Ribas and Stephen 
Wilkins for helpful advice and useful discussions and 
Michael Rauch for inspiring the online tool during 
EWASS 2014. P.B. is supported by European Research 
Council grant StG2010-257080. 
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Fig. 6. — Same as [ 5 ] but here we varied the resolving power R (see Appendix [Aj . Our inferred constraints generally depend weakly on R 
for R IO 1 . For lower resolving power our constraints start to deteriorate faster. 

Zheng, Z., Cen, R., Trac, H., & Miralda-Escude, J. 2011, ApJ, Zheng, Z., & Miralda-Escude, J. 2002, ApJ, 578, 33 
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APPENDIX 

VARYING THE RESOLVING POWER R 

For the main results we presented in the paper we assumed R = 6000. Here, we quantify the impact of varying 
the spectral resolution/resolving power while keeping the SNR fixed at ~ 35. Note that we generate a different noisy 
random realisation of a given shell model for each R. This explains why our constraints do not vary monotonically 
with R for a given model. 
























